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ABSTRACT 

We present a new method for time-efHcient and accurate extraction of the power 
spectrum from future cosmic microwave background (CMB) maps based on properties 
of peaks and troughs of the Gaussian CMB sky. We construct a statistic describing 
their angular clustering - analogously to galaxies, the 2-point angular correlation 
function, ^u{0)- We show that for increasing peak threshold, u, the ^^{0) is strongly 
amplified and becomes measurable for v >1 on angular scales < 10°. Its amplitude 
at every scale depends uniquely on the CMB temperature correlation function, C{9), 
and thus the measured can be uniquely inverted to obtain C{6) and its Legendre 
transform, the power spectrum of the CMB field. Because in this method the CMB 
power spectrum is deduced from high peaks/troughs of the CMB field, the procedure 
takes only [f{u)]'^N'^ operations where /{v) is the fraction of pixels with \5T\ > v 
standard deviations in the map of A'" pixels and is e.g. 0.045 and 0.01 for v='l and 2.5 
respectively. We develop theoretical formalism for the method and show with detailed 
simulations, using MAP mission parameters, that this method allows to determine 
very accurately the CMB power spectrum from the upcoming CMB maps in only 
~ (10"^ - 10^^) X iV^ operations. 

Subject headings: cosmology - cosmic microwave background - methods: numerical 
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1. Introduction. 

By probing the structure of the last scattering surface, the current and upcoming baUoon and 
space borne missions promise to revolutionize our understanding of the early Universe physics. 
This requires probing the angular spectrum of the cosmic microwave background (CMB) with high 
precision on sub-degree scales, or angular wave-numbers I > 200. For cold-dark-matter (CDM) 
models based on inflationary model for the early Universe and adiabatic density perturbations, 
the structure of the CMB should show the signature of acoustic oscillations leading to multiple 
(Doppler) peaks. The relative spacing of the Doppler peaks should then reflect the overall 
geometry of the Universe, whereas the amplitude of the second (and higher) peaks depends 
sensitively on other cosmological parameters, such as the baryon density, re-ionization epoch, etc. 
The recent balloon-borne measurements (de Bernardis et al., 2000, Mauskopf et al 2000, Hanany 
et al., 2000) strongly imply a flat cosmological model because the first Doppler peak occurs at 
I ~ 200 (Kamionkowski, Spergel & Sugiyama 1994, Melchiori et al 2000, Jaffe et al 2001). 

A major challenge to understanding these and future measurements is to find an efficient 
algorithm that can reduce the enormous datasets with A'' ~ 10^ pixels in balloon experiments 
to ~ 3 X 10^ for MAP band with the 0.2° beam (at 90 GHz) to ~ 10^ for the Planck HFI data. 
Traditional methods require inverting the covariance matrix and need ~ operations making 
them impossible for the current generation of computers. Thus alternatives have been developed 
for estimating the CMB multipoles in 0(A^^) operations from Gaussian sky (Tegmark 1997; Oh, 
Spergel & Hinshaw 1999) and general CMB sky (Szapudi et al 2000) as well as study statistics of 
the various methods (Bond et al. 2000; Wandelt et al 2001). 

In this letter we suggest a novel, accurate and time-efficient method for computing the angular 
power spectrum of the CMB temperature field from datasets with large N using the properties of 
the peaks (and troughs) of the CMB field. The peaks are much fewer in number than N, but they 
will be strongly clustered. In Sec. 2 we construct a statistic describing their angular clustering - 
the 2-point angular correlation function ^ in analogy to the galaxy correlation function (Peebles 
1980), i.e. the excess probability of finding two peaks at a given separation angle. We show that 
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this statistic is strongly amplified over the scales of interest (<10°) and should be measurable. 
The value of for a given peak threshold \6T\ > ua would be uniquely related to the correlation 
function of the temperature field C = {6T{x)6T{x+9)). The measurement of ^ can then be uniquely 
inverted to obtain the underlying C and its Fourier transform, the power spectrum Q. This can 
be achieved in just [/(z^)]^A^^ operations, where /(i^) is the fraction of pixels with \6T\ >va and 
is e.g. 4.5-1% for i/=2-2.5. Sec. 3 shows concrete numerical simulations for the MAP 90GHz 
channel in order to estimate cosmic variance, sampling uncertainties, instrumental noise etc. We 
show that with this method the CMB power spectrum is recovered to accuracy comparable to or 
better than by other existing methods but in a significantly smaller number of operations. On an 
UltraSparc II 450 MHz processor the entire sky map with 0.2° angular resolution can be analyzed 
and Cis recovered in only 15 mins and 2.25 hours CPU time for v =2.5 and 2.1 respectively. 

2. Method 

For Gaussian ensemble of N data points (e.g. pixels) describing the CMB data 5 = T — (T) 
one expects to find a fraction f{u) =erfc(i//\/2) with \5\ > ua, where cj^ = (5^) is the variance of 
the field and erfc is the complementary error function. E.g. /(zy)=(4.5, 1, 0.1)xl0-2 for z/=(2,2.5,3) 
respectively. The joint probability density of finding two pixels within d6i^2 of 5i,2 and separated 
by the angular distance 6 is given by the bivariate Gaussian (the vector d has components [6i, 62])- 

11 1 f°° f°° 1 

= ^^^\ '"P^- 2^ • • = (^ loo loo '"P^"^^ ■ '"P^" 2^ ■ ^ ■ ^^^'^ 

where C is the covariance matrix of the temperature field. We model the covariance matrix in (|l|) 
as C{6) = C^dij + C{6ij){l — 5ij), where 5ij is the Kronecker delta and Cq = C(0) +0"^; o"n is the 
noise contribution. We assume that the noise is diagonal, the entire CMB sky is Gaussian and the 
cosmological signal whose power spectrum we seek is contained in C{9) = X](2^ + l)C'/-P/(cos 9) /in. 
The total dispersion of the temperature field is then a = ^C(O) + a'^ 

The peaks of a Gaussian field should be strongly clustered (Rice 1954, Kaiser 1984, Jensen 
& Szalay 1986, Kashlinsky 1987). The angular clustering of such regions can be described by 
the 2-point correlation function, i.e. the excess probability of finding two events at the given 
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separation. The probability of simultaneously finding two temperature excursions with \6T\ > ua 
in small solid angles dwi^2 is dPi2 0c(l + S,)dwidw2- The correlation function of such regions is: 

2CC[p(5l,'^2)+p(-5l,52)]d5ld52 



1 



(2) 



The numerator follows from considering contributions from correlations between 5i , 82 in regions of 
1) 5i > z^v^Co, 82 > i^VCq] 2) 61 < — i/^/Co, 62 < — J^-v/Co; and 3) twice the contribution of Si > i^-v/Cq, 
62 < —i'\/Cq. The "2" in the denominator comes because we consider both peaks and troughs. 

In order to evaluate (|2|) directly, we expand exp[— gig2C'(^)] = J2k^=o ^^'^kP^ (1) 
use the fact that J^^ex.p{—ixy)F{x)x^dx = i^{d^/dy^) J^^exp{—ixy)F{x)dx (Jensen & Szalay 



1987, Kashlinsky 1991). Because J^^exp{-iq5) exp{-q'^Co/2)dq = y/27r/Co exp{-^) we get: 
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Substituting ^ into (|2|) allows to expand ^,^{0^ into the Hermite polynomials, 
Hn{x) = {—y- exp{x'^){d"- /dx^') exp{—x'^) , to obtain: 

C . 



with: 
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where H-i{x) = ^ exp(x^)erfc(x). At each angular scale the value of for every v is determined 
uniquely by C at the same 0. Note that in the limit of the entire map (i^=0) our statistic (,u=0 
and our method becomes meaningless; the new statistic has meaning only for sufficiently high u. 
One should distinguish between the 2-point correlation function, ^, we directly determine from the 
maps, and the commonly used statistics in CMB studies, the temperature correlation function, C. 

Fig. 1 shows the properties of the left panel shows the variation of with C/Cq for 
fixed and the middle panel shows the variation of (^i, with for fixed C/Cq. The first term 
in the sum in eq. (^) contains H'Kocu'^), but as the middle panel of Fig.l shows for sub-degree 
scales (where | C/Cq | > 0.1) changes more steeply than u'^. This means that k>l terms are 
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important for accurate inversion of in terms of Q's. The right panel shows vs the angular 
separation 9 for v=2 for two flat CDM models: ACDM model (thin line) with (fitotab ^a)=(1)0.7) 
and Obaryon^^=0.03 and SCDM with (Ototab f^A)=(l,0) and r2baryon/i^=0.01 (thick line). The 
first model has prominent Doppler peaJcs and baryon abundance in agreement with BBNS, the 
second model requires significantly higher baryon abundance but has a much smaller second 
Doppler peak. There would be non-linear to quasi-linear (and easily detectable) clustering of high 
peaks out to the angular scale where C{6) drops to only ~0.1 of its maximal value at zero-lag. 
This covers the angular scales of interest for determining the sub-horizon structure at the last 

— 1/2 

scattering. Because the uncertainty in measuring ^ is ~ 

-^pairs (Peebles 1980), the value of ^ can 
be determined quite accurately in non-linear to quasi-linear regime. At the same time, as the left 
panel in Fig.l shows, over this range of scales the amplitude of ^i, changes rapidly with C making 
possible a stable inversion procedure to obtain C{6) from 

This suggests the following procedure to determine the power spectrum of CMB in only 
c^p{y)N'^ operations: • Determine the variance of the CMB temperature, Co, from the data in 
N operations; • Choose sufficiently high v when f{v) is small but at the same time enough pixels 
are left in the map for robust measurement of ^jy(0); • Determine E,u{&) in [/('^)]^-^^ operations. • 
Finally, given the values of {Cq,v) solve equation Ai,{C / Cq)=^v{S) to obtain C{0) and from it C;. 

3. Numerical results and applications 

In order to apply the proposed method in practice and to estimate the cosmic variance, 
sampling and other uncertainties, we ran numerical simulations with parameters corresponding 
to the MAP 90 GHz channel. The CMB sky was simulated using HEALPix (Gorski et al 1998) 
software with Afsidc=512 and Gaussian beam with FWHM=0.21° (or /Nyqmst=640) for the SCDM 
and ACDM models. To this we added Gaussian white noise with the rms of 35/xK per 0.3° x 0.3° 
pixel (Hinshaw 2000). We assumed that the foreground contribution at 90 GHz can be subtracted 
to within a negligible term. Fig. 2 shows the SCDM model sky with both peaks and troughs with 
\5\ > u\/Cq, v=2, marked with white dots. The clustering of peaks/troughs is very prominent 
especially on small scales, but the clustering pattern is very different between the models. 
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To determine we divided 9 into 31415 equally spaced bins between and 180° and the 
number of pixel pairs, A^i2, with \5\ > va, was oversampled and determined in each bin. This is 
the dominant CPU time-consuming procedure of {f{v)N]^ operations. The " raw" value of ^raw 
was determined in each of the 31415 bins as Craw=-^i2/-^rr — 1; where N^-r is the number of pairs 
for a Poissonian catalog with the total number of pixels equal to the number of peaks and troughs 
in our CMB sky. The final ^ was determined as follows: the angular interval between and 180° 
was divided into bins centered on the roots of the 800th order Legendre polynomial in order to 
facilitate the later inversion of C{9) into C/'s via the Gauss-Legendre integration. The final ^ was 
obtained from ^raw by convolving the latter with a Gaussian filter of 4' dispersion centered on each 
of the 800 Legendre polynomial roots. The value of thus obtained ^(0) is shown for one realization 
of the two CDM models in the right panel of Fig.l; it agrees well with eqs. (^^). 

Having fixed v and determined Co from the map we now solve the equation ^y{6)=A^{C {9) / Cq) 
with respect to C{9) with Ay given by eq.(^ and determine C{9) at each of the roots of /=800 
Legendre polynomial. In the final step the multipoles were determined by direct Gauss-Legendre 
integration of Ci=2tt J C (9) Pi{cos 9) sin 9 d9. At 9 >10°, where is very small and hard to 
determine (see Fig.lc), our recovered C{9) has larger uncertainties. Because we are interested in 
high I multipoles the recovered correlation function C{9) was further tapered above 15°. (At / of 
interest the results are insensitive to details of tapering). To check the statistical uncertainties in 
the determination of Q's we ran 700 simulations for z^=2.5 and 350 simulations for 1^=2.1. 

Fig. 3 shows the results of the numerical simulations for SCDM model with the instrument 
noise of the MAP 90 GHz channel. The distribution of the multipoles determined from the 
simulated maps in this method is shown in Fig. 3a, b for /=200 (the first Doppler peak), /=350 
(the first trough) and /=475 (the second Doppler peak for SCDM model). The best-fit Gaussians 
are shown with smooth lines; they give good fit to the histograms at high I. The 68% confidence 
limits on Q's are very close to the dispersion of the best-fit Gaussians shown and the 95% limits 
are roughly twice as wide for all Vs as would be the case for approximately Gaussian distributions. 

Fig. 3c shows C{9) determined by our method from one realization for i/=2.1 (solid line) 
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and i/=2.5 (dashes). Dotted line shows the theoretical C(0)=X^(2Z + l)QP/(cos^)/47r with 
SCDM values of Cfs. The C{0) determined with our method is within 5% of the theoretical 
value. This uncertainty is within the cosmic variance of small-scale C{9) from low-/ contribution 
(predominantly quadrupole) which is ~ 500/LtK2 (Bennett et al 1996, Hinshaw et al 1996). 

Finally, Fig. 3d shows our full sky determination of the CMB power spectrum from the 
synthesized SCDM maps. Solid line shows the theoretical (input) spectrum juxtaposed with the 
power spectrum determined with the peaks method for z/=2.1 (diamonds) and 2.5 (crosses). The 
latter was band-averaged into AZ=50 wide bins and the symbols are plotted at the central bin value 
+/- 5 for i/=2. 5/2.1 respectively to enable clearer display. Band power averaged over larger Al will 
have less variance but will give fewer independent data points. In computing the band-averages 
we gave equal weight to all multipoles. This top-hat window function is not an optimal estimator 
of the power spectrum (Knox 1999). However, we checked that with this window function, the 
multipoles at the different Z-bins were not correlated. The shown uncertainties correspond to the 
dispersion from the Gaussian fits to the distributions and are very close to the 68% error bars 
which are approximately half the 95% errors. With our method for u=2.1 we recover the power 
spectrum with variance only 1% larger than the full-sky cosmic variance limit at the central I at 
the first Doppler peak (Z=200), 2% larger at /=350 and 8% larger at the second peak (/=475). We 
compare the uncertainty of our numbers to the full-sky cosmic variance which is, excluding the 
noise, ACl=^/2/{2^+TjCl. Note that for our method we simulated the two-year full-sky MAP 90 
GHz channel, so our assumed noise is ~100/iK per 7' pixel with the FWHM=12.6' beam, i.e. twice 
that assumed in Szapudi et al (2000) with the FWHM=10' beam who obtain a similar accuracy 
with direct computation of C{0) for a small BOOMERANG-size patch of the sky. At Z<ZNyquist of 
the beam our method determines Q's without bias. The accuracy of our method for a given Al 
can be further improved by going to ^'<2.1, but increasing somewhat computational time (e.g. for 
i/=2 the computational time will increase by 60% over v=2.1). 

4. Conclusions 

We introduced here a new statistic, ^jy, and showed that with it one can recover the power 
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spectrum from Gaussian CMB maps in a very accurate and time-efficient way. We have shown 
that for peaks and troughs of such temperature field, their angular 2-point correlation function, 
^i/, would be strongly amplified with increasing threshold v and would be measurable. Because 
its amplitude at a given angular scale depends on the amplitude of the temperature correlation 
function, C, at the same scale, the former is then inverted to obtain the power spectrum of the 
CMB. The method requires [f{i')N'^<^N'^ operations. For balloon experiments it would work for 
z/~l-1.5 or in (10~^ — 10~^)iV^ operations, for MAP highest resolution with ~3xl0^ pixels it can 
work at z^~2-2.5 or in (2xl0~^ — 1.5 x 10^^)A^^ operations and for higher resolution maps, such as 
Planck HFI maps, still higher u can be used leading to accurate results in < 10~^A/^^ operations. 
We demonstrated with simulations that for the two-year noise levels for MAP 90GHz channel with 
this method we can recover the CMB multipoles with z/=2.1, or in 1.2xlO~^iV^ operations, out to 
^Nyquist with Uncertainty only a few percent larger than the full-sky cosmic variance. We assumed 
a diagonal noise covariance matrix, but the method can be extended to another Gaussian noise 
provided its covariance matrix is known. Because here we work with the correlation functions, our 
method is immune to geometrical masking effects e.g. from Galactic cut and other holes in the 
maps. This means that we can remove regions where the foregrounds are very bright, or the noise 
inhomogeneity is too high. With the exception of isolated areas, the MAP noise variations are 
expected to be ~10-15% over the cut sky (Hinshaw, private communication). Because cr^<CC(0) 
such variations would lead to only small variations in the effective v. The method applies to 
Gaussian CMB sky which is expected in conventional models. If the CMB sky turns out to be 
non-Gaussian, our method may not be applicable, or may need substantial modifications by e.g. 
looking at for peaks and troughs separately. 

We acknowledge fruitful conversations with Gary Hinshaw and thank Keith Feggans at NASA 
GSFC for generous computer advice and resources. C.H.M. and F.A.B. acknowledge support of 
Junta de Castilla y Leon (project SA 19/OOB) and Ministerio de Educacion y Cultura (project 
BFM2000-1322). 
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FIGURE CAPTIONS 

Fig. 1: (a) vs C/Cq for u = 0.5,1,1.5,2,2.5,3 from bottom to top. (b) ^i, vs v for 
C/Cq = 0.1,0.2,0.3,0.5,0.75,0.95 from bottom to top. (c) vs 9 for u = 2 m. one realization of 
the two CDM models: Plus signs correspond to ^ determined directly from one simulated map of 
ACDM with FWHM=0.21° resolution and the noise corresponding to the 90 GHz MAP channel; 
diamonds show the same for SCDM. Thick and thin solid lines show the values of ^{9) from eqs. 
(HI) for SCDM and ACDM respectively. 

Fig. 2: All sky distribution of pixels with v=2 for SCDM model. 

Fig. 3: (a), (b) Histograms of the recovered C;'s for Z=200, 350 and 475 are shown for 
1^=2.5 (top) and 2.1 (bottom). Smooth lines show the best-fit Gaussians to the histogram data, 
(c) C{0) vs 6 for SCDM model: theoretical value is shown with dotted line. The values for one 
realization are shown with solid (for u=2.1) and dashed (z^=2.5) lines, (d) C; vs / for SCDM 
model. Solid line corresponds to the theoretical input value. The spectrum recovered with our 
method from simulated 90 GHz MAP maps is shown after band-averaging with AZ=50 with filled 
diamonds {u=2.1) and crosses (i^=2.5). To enable a clearer display the central values of multipoles 
are shifted by 5 to the left for v=2.\ and to the right for u=2.\. The error bars correspond to 
the dispersion of the Gaussian fits such as as shown in Fig. 3a and practically coincide with 68% 
confidence limits which in turn are approximately half the 95% limits. 
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Fig. 2.— 



